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ABSTRACT 

Aims. We examine radial and vertical metallicity gradients using a suite of disk galaxy hydrodynamical simulations, supplemented 
with two classic chemical evolution approaches. We determine the rate of change of gradient slope and reconcile the differences 
existing between extant models and observations within the canonical "inside-out" disk growth paradigm. 

Methods. A suite of 25 cosmological disks is used to examine the evolution of metallicity gradients; this consists of 19 galaxies 
selected from the RaDES (Ramses Disk Environment Study) sample (Few et al., in prep), realised with the adaptive mesh refinement 
code ramses, including eight drawn from the 'field' and six from 'loose group' environments. Four disks are selected from the 
MUGS (McMaster Unbiased Galaxy Simulations) sample (Stinson et al. 2010), generated with the smoothed particle hydrodynamics 
(SPH) code gasoline, alongside disks from Rahimi et al. (201 1 : gcd+) and Kobayashi & Nakasato (201 1 : grape-SPH). Two chemical 
evolution models of inside-out disk growth (Chiappini et al. 2001; Molla & Diaz 2005) were employed to contrast the temporal 
evolution of their radial gradients with those of the simulations. 

Results. We first show that generically flatter gradients are observed at redshift zero when comparing older stars with those forming 
today, consistent with expectations of kinematically hot simulations, but counter to that observed in the Milky Way. The vertical 
abundance gradients at ~ 1-3 disk scalelengths are comparable to those observed in the thick disk of the Milky Way, but significantly 
shallower than those seen in the thin disk. Most importantly, we find that systematic differences exist between the predicted evolution 
of radial abundance gradients in the RaDES and chemical evolution models, compared with the MUGS sample; specifically, the 
MUGS simulations are systematically steeper at high-redshift, and present much more rapid evolution in their gradients. 
Conclusions. We find that the majority of the models predict radial gradients today which are consistent with those observed in 
late-type disks, but they evolve to this self-similarity in different fashions, despite each adhering to classical 'inside-out' growth. We 
find that radial dependence of the efficiency with which stars form as a function of time drives the differences seen in the gradients; 
systematic differences in the sub-grid physics between the various codes are responsible for setting these gradients. Recent, albeit 
limited, data at redshift z~1.5 are consistent with the steeper gradients seen in our SPH sample, suggesting a modest revision of the 
classical chemical evolution models may be required. 
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1. Introduction 

The recognition that metals are not distrib uted homogeneously 
throughout the disk of the Milky Way (IShaver et all [l983) 
has proven to be fundamental in our efforts to understand 
the role of interactions, mergers, accretion, migration, and 
gas flows, in shaping the formation and evolution of galax- 



These authors contributed equally to this work. 



ies. A rich literature now exists which confirms these radial 
abundance trends in both spirals (e.g 
Afflerbach et all [l997; Molla et al 



Simpson et al. 1995 
1999t ICarreraetal l 2008 



rtinciuacii ci ai. iyy i , iviuna ci ai. Lyyy, v^aucia ci ai. zuuo, 

iKewlev et alJl2010t ISanchez-Blazquez et al]' 201 lh and e llipti- 



cals ( e.g. iKormendv & D iorgovskil ll989t iFranx & I llingworth 
1990; iPeletier et al] 119901) . Verti cal trends have been stud- 
ied somewhat less frequent l y (e.g. iMarsakov & Borkoval '2005. 
2006; ISoubiran et al] 12008; Nav arro et al] 1201 ll) . but provide 
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unique insights into the discrete nature (or lack thereof) of the 
thin disk - thick disk interface (and associated kinematical heat- 
ing processes). 

Observations of nearby spiral galaxies show that the inner 
disks have higher metallicities than their associated outer disk 
regions; at the present day, typical gradients of ~-0.05 dex/kpc 
are encountered. These somewhat shallow gradients have pro- 
vided critical constraints on models of galaxy formation and 
evolution, and are fundamental to the predictions of the classi- 
cal "inside-out" paradigm for disk growth. Predictions have been 
made of the time evolut ion of metallicity gradients in c hemical 
evolution models (e.g. lMolla et alJfl997blFu etal. 2009) an d ob - 
servationally from plantetary nebulae (e.g. Maciel et al. 2003), 
although until recently, we have had essentially no direct obser- 
vational constraints on what the magnitude of the time evolu- 
ti on of the gradients should be. This has c hanged with the wor k 
of lCresci et all d2010b Hones et alj fcOlOh . lOuevrel et alj ( 1201 lb . 
and Yuan et all ( 201 lb . who have, for the first tim e, extended 
radial abundance gradient work to high redshifts. lYuan et al] 
(1201 lh show that for at least one "Grand Design" disk at 
redshift z~1.5, the metallicity gradient is significantly steeper 
(-0.16 dex/kpc) than the typical gradient encountered todayQ 
Constraining the metallicity gradients of galaxies beyond the lo- 
cal Universe remains a challenge for the future. 

Usi ng SPH simulations of disk galaxy mergers, Rupke et al. 
(2010a|) show strong correlations of metallicity with environ- 
ment and merger history, focussing on the effec ts of gas in- 
flows and star formation rate. Observations by ICooper et alJ 
(2008) show that higher metallicity galaxies are more abundant 
in group enviroments and lKewlev et alj d2006b showed that inter- 
acting pairs of galaxies have systematically lower metallicities 
(~0.2 dex lower) than field galaxies or more loosely associated 
pairs. Radial gradients have been sho wn to flatten for g alaxies 
that have experienced recent mergers dKewlev et alJ2010l) : these 
also result in higher velocity disper sions and redistributio n of the 
cold gas. In agreement with this, Mic hel-Dansac et ail (2008) 
studied the mass-metallicity relation for merging galaxies and 
concluded that the infall of metal poor gas during merger events 
lowers the gas phase metallicity. However, the timescale over 
which redistributed gas develops into a gradient like those we 
see in spiral galaxies today is unknown. 

There have been several studies of chemistry within cos- 



mological hydrodynamical 
Kawata & Gibson 



1996 



simulations (e.g 



2003 



Raiteri et al. 
Okamoto etal] 12008 



Scannaoieco et alJ " 2008t IZolotov et al.1 201 J~ Rahimi et alJ 
2010t IWiersma et alj|201 it iKobavashi & Nakasatoll2oT 1), each 
modelling certain observational properties with varying degrees 
of success. Some studies have examined the radial and/or 



vertica l gradients using hyd rodynamical codes (e.g. Rupke et al. 
1201 Oat iRahimi et al.1 120 l ib , but the numerical study of radial 
gradients has predominantly been in the context of classical 
galac ti c chemical evo l ution c odes (e.g. iPrantzos & Boissierl 
2000; IChiappini et al] 1200 It iMolla & Dfazl 12005b . In this 
paper, we use 25 simulations realised with three different 
cosmological hydr odynamical codes: gaso line dWadslev et alj 
2004) and gcd+ dKawata & Gibsonl l2003h . both g ravitational 
N-Bo dy + Smoothed Partic le Hydrodynam ic (SPH) (Monaghan 
11992b codes, and ramses (Teyssier 2002), an Adaptive Mesh 
Refinment (AMR) code. Alongside th ese, we use the results 
from the chemical evolu tion models of Chiap pini et alj d2001b 
and lMolla & Dfazl d2005b . 

Our work aims to fill an important gap i n the field, by com- 
plementing orbita l parameter studies (e.g. Ru pke et alj |2010a; 
IPerez i 



nting orbital parameter studies (e.g. Kupke et al. 2UlUa; 
et alJ 2011b. systematic sub-grid physics parameter stud- 



ies (e.g. Wiersma et alj 20ll|). and detailed dissections of in- 
dividual systems (e .g. Rahimi et ai] |201 It IZolotov et al.1 1201 Ot 
IKobavashi & N akasato 201 1), with a statistical sample of Milky 
Way-like analogs. Our approach is different, but c omplementary, 
to the careful and compelling parameter study of IWiersma et al.1 
d2011b : their, the goal was to vary the input physics and examine 
the outcome, regardless of whether or not the simulated end- 
products might be classified still as Milky Way-like. Instead, 
we have sampled a range of codes, sub-grid physics, and initial 
conditions, each of which has been 'calibrated', in some sense, 
by their respective authors, to resemble a classical Milky Way- 
like system. With that calibrated sample, our unique contribution 
is to examine the 'path' by which the gradients evolve, search 
for both random and systematic trends/differences between the 
samples, and compare with new empirical data at high-redshiftQ 
This is the first time such a comparison of the temporal evolution 
of metallicity gradients has been undertaken with a statistical 
sample of simulated disk galaxies. 

The outline of the paper is as follows. The main differences 
between the codes are described in ^2j where we concentrate pri- 
marily upon the relevant mechanisms associated with the treat- 
ment of star formation and feedback (both energetic and chemi- 
cal). The metallicity gradients inferred today for stellar popula- 
tions of different ages are presented in $3] This is expanded upon 
in $4] where the radial metallicity gradients of the young stellar 
population as a function of redshift are considered. Finally, we 
summarise our findings in §0 



2. Simulations 



1 At even higher redshifts (z~3.3), Cresci et al. (2010) and Troncoso 
et al. (2012, in prep), as part of the AMAZE/LSD surveys, suggest 
that both inverted gradients (higher abundances in the outskirts, rel- 
ative to the inner disk) and standard declining gradients are seen. 
From the latter surveys, inverted gradients (ranging from +0.0 to 
+0.1 dex/kpc) appear associated with very massive stellar disks at 
these high-redshifts (M„>3xl0 9 M Q ), while declining gradients (rang- 
ing from -0.0 to -0.2 dex/ kpc) appear associa ted with lower mass stel- 
lar disks (M,<3xl0 9 M^I. ICresci et al.1 d2010l) suggest that the inverted 
gradients are due perhaps to recent infall of pristine material into the 
inner disk. These Lyman Break Galaxies, with their ~l-2 orders of 
magnitude greater star formation rates (relative to the typical Milky 
Way progenitor at that redshift), are mo re likely associa ted with massive 
spheroids in clusters/groups today (e.g. Naeamine 2002), as opposed to 
the Milky Way, and so are not directly comparable with the simulations 
described here. 



The simulations used in this paper are fully described in Stinson 
et al. (2010: MUGS), Rahimi et al. (2011: Gall), Kobayashi & 
Nakasato (2011: KN11) and Few et al. (in prep: RaDES); the 
main characteristics of the simulations and their parent codes 
are described here and itemised in Table [T] The chemical evo- 
lution models are ful ly described in Chiap pini et alJ d200lb and 
IMolla & Dfazl d2005b . but again we describe the main aspects in 
the following section. 



2 In spirit, this is exactly the approach taken in th e seminal Galactic 
Chemical Evolution Comparison Project (Tosi 1996), which examined 
the time evolution of classic chemical evolution models calibrated to 
the solar neighbourhood, in order to see where they differed 'away' 
from this calibrated boundary condtion. 
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Table 1. Basic present-day (z=0) characteristics of the 25 simulated disks. Column (1): simulation suite to which the the code used 
to simulate the galaxy (Column (2)) belongs; Column (3): total (dynamical) mass within the virial radius; Column (4): mass of the 
stellar disk, after application of the kinematic and spatial cuts described in § 3; Column (5): exponential scalelength of the stellar 
disk; Column (6): local environment of the galaxy; Column (7): mass-weighted vertical stellar abundance gradient, averaged over 
the radial range 1.5<r ( ii s i c 2.5; Column (8): mass-weighted radial young (stars born within the past 100 Myrs) stellar abundance 
gradient, after application of the kinematic and spatial cuts described in § 3. 



2.1. MUGS 

The MUGS galaxies were run using the gravitational N-body 
+ SPH code gasoli ne which was introduced and described in 
Wads lev et alJ (2004). Below, we emphasise the the main points 
concerning the star formation and feedback sub-grid physics 
used to generate this suite of simulations, but first remind the 
reader of the background framework in which they were evolved, 
in addition to their basic characteristics. 

The MUGS sample (IStinson et all |2010) consists of 
16 galaxies randomly drawn from a cosmological volume 
50/i _I Mpc on a side, evolved in a Wilkinson Microwave 
Anisotropy Probe Three (WMAP3) ACDM cosmology with 
H = 73 km s" 1 Mpc" 1 , Q„, = 0.24, Q A = 0.76, Q. b = 0.04, and 
cr 8 = 0.76. Each galaxy is resimulated at hi gh resolution by us - 
ing the volume renormalization technique (Klv pin et al.ll200lh . 
with a gravitational softening length of 310 pc. The galaxies 
range in mass from 5xl0 n M to 2xl0 12 M Q . The four galax- 
ies with the most prominent disk^3 were selected: g42zQ, gl536, 
g24334, and gl5784, the latter of which is the closest to a Milky 
Way analog in the sample. 



3 By 'prominent', we mean the inclusion of those for which there 
was unequivocal identification of the disk (from angular momentum 
arguments constructed from the gas and young star distributions, as dis- 
cussed in §3.1. In a secondary sense, this eliminated extreme values 
of bulge-to-total, but formally, we only included those disks for which 
alignment based upon the gas/young stars was obvious. 

4 g 422 was not described in the original MUGS paper (Stinso net alj 
2010); it was produced identically to the MUGS suite and will be de- 
scribed fully in an upcoming paper. 



Sta r formation and sup ernovae feedback uses the blastwave 
model ( St inson et al.ll2.006h whereby gas particles can form stars 
when they are sufficiently dense (>1 cm 4 ) and cool (< 15000 K). 
Gas particles which satisfy these criteria can form stars accord- 
ing to the equation ^jf--c* ■^ 2i , where c* is the star formation 
efficiency and is fixed to be 0.05. M gas is the mass of the gas par- 
ticle forming the star particle of mass M+ and td yn is the dynam- 
ical time of the gas. Heatin g from a uniform ultra violet ionising 
background radiation field (Haardt & Mada ui 19961) is employed, 
and cooling is derived from the contributions of both primordial 
gas and m etals; the metal cool ing grid is derived using CLOUDY 
(v.07.02: lFerland et al.l (Tl998)). u nder the as sumption of ionisa- 
tion of equilibrium, as detailed by Shen et al. (201§. 

The c hemical evolution m odel used in gasoline is fully de- 
scribed in Rait eri et al.l (Il996al) : here, we only discuss the main 
points. All stars with masses above 8 M explode as Type II su- 
pernova (SNell). An efficiency factor couples 40% of a given 
supernova's energy (10 51 erg) to the surrounding interstellar 
medium (ISM). The metals that are tracked in this version of 
gasoline (O and Fe) all come from supernovae and a re allowed 
to diff use between neighbouring SPH particles, after lShen et al.l 
d2010b . The Type la supernovae (SNela) eject iron and oxygen; 
for every SNIa, 0.76 M of 'metals' is ejected, divided between 
iron (0.63 M Q ) and oxygen (0.13 M ). Our binary model for 
Type la supernovae is based upon the s ingle- degenerate progen- 
itor formalism of Greg gio & Renzinil ( 19831) . with secondaries 
spanning in mass from 1.5 to 8.0 M |^ Enrichment from SNell 



5 We have excluded secondaries in the 0.8 - 1.5 M G range; doing so, 
regardless of IMF, only impacts on the SNela rate at the ~20% level. 
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is based upon po wer law fits in stellar mass t o the nucleosynthe- 
sis yield tables of | Wooslev & Weaver! (1 19951) . convolved with a 
Kroupa dKroupa et al.ll 1993b initial mass function (IMF), in or- 
der to determine the mass fraction of metals ejected. The total 
metallicity in this version of the code is tracked by assuming 
ZsO+FeQ For these runs, only the Z=Z yields were used, and 
long-lived SNela progenitors (those with secondaries with mass 
m< 1 .5 M ) were neglected. 



2.2. Gall 

Gall is a h igher-resolution re-simulation of galaxy Dl 
from iKawata et all (|2 004) using the SPH code gcd+ 
dKawa ta & Gi bsonll2003l) : while its charac t eristics have 
been d iscuss ed previously b y iBailin et al.l (2005), Rahimi et al.l 
d2010l) . and iRahimi et al.l (1201 ll) . an overview is provided 
here for completeness. Employing a comparable volume 
renormalisation / 'zoom-style' technique to that described in 
§ 2.1 (with a gravitational softening of 570 pc in the highest 
resolution region), Gall was realised within a ACDM cosmo- 
logical framework with Hq = 70 km s _I Mpc~', £2 m = 0.3, 
Q A = 0.7, Qi = 0.04, and cr 8 = 0.9, resulting in a Milky 
Way analog of virial mass 8.8xl0 n M . The effect of the 
ultraviolet background radiation field was neglected, while 
metal -dependent radiative coo ling (adopted from MAPPINGS- 
III dSutherland & Dopitall 19931) ) was included. 

The star formation prescription employed requires (i) the gas 
density to be above a threshold of 0.1 cm 4 , (ii) a convergent 
gas flow to exist, and (iii) the gas to be locally Jeans unstable. 
A standard Salpeter initial mass function (IMF) was assumed, 
along with pure thermal feedback from both SNell and SNela 
(10 50 erg/SN) being coupled to the surrounding SPH particles. 

The chemical evolution implementation within GCD+ 
takes into account the metal-dependent nucle osynthetic 
byproducts of SNell dWooslev & Weaverl 1 19951) . SNela 
(Iwam otoetalJf l999). and low- and intermediate-mass AGB 
stars dvan den Hoek & Groenewegenl Il997l) . Relaxing the 
instantaneous recycling approximation, GCD+ tracks the 
temporal evolution of the nine dominant isotopes of H, He, C, 
N, O, Ne, Mg, Si, and Fe. The SNela progenitor formalism of 
iKobavashi et all d2000h is adopted. 



2.3. KN1 1 

KN11 corr esponds to the so-called 'Wid er Region' model de- 
scribed by IKobavashi & Nakasatoj d201 lb . realized used a hy- 
brid grape- SPH code. This model was drawn from the 5 Milky 
Way-analogs which eventuated from a larger suite of 150 semi- 
cosmologicaQ simulations. The cosmological parameters em- 
ployed match those of §2.2, and led to a Milky Way analog 
of mass l.lxlO 12 M . The effect of the ultraviolet background 
radiation field was included, as was metal-dependent radiative 



6 By assuming Z=0+Fe, we admittedly underestimate the global 
metal production rate by nearly a factor of two; our next generation runs 
with gasoline employ a more detailed chemical evolution model, in- 
corporating the nucleosynthetic byproducts of asymptotic giant branch 
evolution and thereby ameliorating this effect. 

7 By 'semi-cosmological', we mean that the simulated field was 
not large enough to sample the longest waves (and, as such, underes- 
timate the degree of gravitational tidal torque which would otherwise 
be present in a fully cosmological framework), and so the initial sys- 
tem is provided with an initial angular momentum via the application 
of rigid rotation with a constant spin parameter A=l. 



coolin g (adopted from MAPPINGS-III dSutherland & D onita 
119931) 1. 

The star formation prescription employed requires (i) the 
gas density to be cooling, (ii) a convergent gas flow to exist, 
and (iii) the gas to be locally Jeans unstable. The star formation 
timescale is chosen to be proportional to the dynamical timescale 
(f sf =f dyn /c), where the star formation efficiency is chosen to be 
c=0.1. A standard Salpeter initial mass function (IMF) was as- 
sumed (with lower and upper mass limits of 0.07 and 120 M , re- 
spectively), along with pure thermal feedback from both SNelQ 
and SNela (~10 51 erg/SN) being distributed to the surrounding 
SPH particles within 1 kpc (weighted by the SPH kernel). 

The chemical evolution implementation within grape-SPH 
takes into accou nt the metal-dependent nucleosyn thetic byprod- 
ucts of SNell dKobavashi et al.l 120061) . SNela dNomoto et al. 
1997), and low- and intermediate-mass AGB stars (Karakas 



2010). 



2.4. RaDES 

The third galaxy sample (RaDES: Ramses Disk Environment 
Study) was simulate d using the ad aptive mesh refinement 
(AMR) code ramses dTevssier| 120021) . The motivation behind 
these simulations was to determine the systematic differences 
between simulated galaxies with neighbouring dark matter 
haloes similar to the Local Group and those in the field. 
The ramses simulations include gravity, radiative cooling, and 
heating from a uniform ionising UV background radiation 
dHaardt & Mad au 1996). Hydrodynamic behaviour of the gas 
phase and gravitational potential is calculated on a spatially 
adaptive grid. A full de scription of the star format ion model used 
in ramses is given by iDubois & Tevss ier (2008); here we give 
just a brief account of its implementation. 

Gas cells with density greater than a given threshold allow 
stars to form at a rate proportional to the density, p = -p/t+, 
where f* is the star formation timescale, which itself is propor- 
tional to the dynamical t ime (fo(p /po)~ 1 '' 2 ), as first d escribed by 
iRasera & Tevssierl d2006l) . After IDubois & Tevssierl d2008l) . we 
use a threshold of po =0.1 cirT 3 and to = 8 Gyr. In combi- 
nation, these choices correspond to an adopted star formation 
efficiency of 2%. Feedback from SNelQ occurs instantaneously 
and the mass carried away is parameterised as (rjs n + t]w), where 
7/sjv is the fraction of a stellar particle's mass that is ejected by 
SNell and rjw is the fraction that is swept up in the SNII wind. 
In the RaDES simulations, jjsn - 0.1 and rjw - 0, which for 
these runs, led to less strongly peaked rotation curves. Energy is 
injected into the gas phase in the form of thermal and kinetic en- 
ergy, distributed across a superbubble of radius rabble according 
to a Sedov blastwave formalism. The metallicity of SN ejecta is 
determined by converting a fixed fraction, fz, of the non-metal 
content of new stars into metals; all galaxies in the RaDES sam- 
ple used /z=0.1. 

RaDES is comprised of two subsamples allowing for a sta- 
tistical intercomparison of field galaxies and those in environ- 
ments similar to those of loose groups; the full details will be 
presented in Few et al. (in prep). These simulations take place 



8 50% of the massive stars are assumed to end their lives as SNell, 
while the remaining 50% are assumed to end their lives as lOx more 
energetic hypernovae. 

9 SNIa are not accounted for in RaDES, although we have recently 
completed a chemical evolution upgrade to RAMSES which parallels 
that implemented within GCD+ (§ 2.2); this will described in a future 
work. 
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in either 20/i~' Mpc (grid resolution of 440 pc) or 24/T 1 Mpc 
(grid resolution of 520 pc) volumes with 512 3 dark matter 
particles in the central region. The cosmology of these boxes 
is H =70 km s^Mpc" 1 , Q A =0.72, Q m =0.28, Q fe =0.045, and 
o- 8 =0.8. 

The sample employed here consists of nine isolated (field) 
galaxies and ten situated within loose groups. The latter are de- 
fined as being those for which two L, halos of comparable mass 
reside within 1.5 Mpc of one another, and neither are located 
within 5 Mpc of a halo with mass in excess of 5xl0 12 M Q . The 
latter criterion avoids the proximity to rich clusters. In a sta- 
tistical sense, these 'loose groups' can be thought of as Local 
Group analogs, at least in terms of dynamical mass, proxim- 
ity to companion galaxies, and the avoidance of rich clusters. 
The field sample contains those halos that are even more iso- 
lated from neighbouring massive halos : specifically, no Mhaio > 
3xlO n M within 3 Mpc). The virial mass range of the RaDES 
sample spans 2.5x10" to 1.6xl0 12 M Q . 

2.5. Chemical Evolution Models 

In this work, we compare our results from the hydrodynami- 
cal simulations described in § 2.1-2.3 to two chemical evolu- 
tion models both designed to reprod uce the main feat ures of our 
Galaxy. The models a re described by Chiap pini et al.l d2001l) and 
Moll a & Dfazl (120051) . and we refer the reader to these papers for 
full details. 

In the model by IChiappini et alj (1200 ll) . the Milky Way 
forms by means of two main infall episodes, both represented 
by exponential infall rates. The first infall episode, characterised 
by the rate cr H ocA e~^ T "^ B , is associated with the formation of 
the halo and thick disk, with an e-folding timescale (t^^h) of 
~1 Gyr. The constant A is determined by requiring that the 
present-day mass surface density of the halo is reproduced. 

The second infall phase is represented as cfo^BiR) e~'/ T »>/», 
and is associated with the formation of the thin disk. The thin 
disk is represented by independent annuli, each 2 kpc wide, with 
no exchange of matter between them (i.e., no radial gas flows). 
The e-folding timescale (t,„/ ; b) of the second infall is assumed 
to be a linear function with increasing galactocentric radius (i.e., 
Tinf,D(R) re R) - enforcing the so-called "inside-out" paradigm 
for disk growth, with the gas accumulating faster in the inner re- 
gions of the disk, relative to the outer disk. The timescales here 
vary from ~2 Gyr in the inner disk, to ~1 Gyr in the solar neigh- 
bourhood, and up to ~20 Gyr in the outermost parts of the disk. 
The constant B(R) is fixed in order to reproduce the present-day 
total surface mass density (stars + gas) in the solar neighbour- 
hood. The star formation rate cr„ is expressed by the common 
Schmidt-Kennicutt law, cr, <x vcr k gas (R, t), where cr gas (R, t) rep- 
resents the gas density at the radius R and at the time f, and 
k = 1.5. The star formation efficiency v is set to 1 Gyr -1 , and 
becomes zero when the gas surface density drops below a cer- 
tain critical threshold, adopted here to be o~ t h=7 M pc~ 2 . The 
nucleosynthesis prescriptions for AGB stars and SNela+SNell 
are drawn from the same sources liste d in § 2.2. 

The chemical evolution mode l of iMolla & Dfazl (120051) dif- 
fers from that of Chiap pini et al.l d2001l) in several aspects, in 
that it is multiphase, treating the ISM as a mixture of hot dif- 
fuse gas and cold molecular clouds. Each galaxy is assumed to 
be a two-zone system, comprised by a halo formed in an early 
gas-rich phase and a disk. The gas of the disk is acquired from 
the halo through an imposed infall prescription characterised by 
the inverse of the collapse time, which itself depends upon the 
total mass of the galaxy. The mass profile is imposed to ad- 



here to the Per sic e aTHl996) universal rotation curve. Similar 
to IChiappini et alJ ( 20011) . each galaxy is divided into concen- 
tric cylindrical zones 1 kpc wide. The collapse timescale de- 
pends on radius via an exponential function t( R) qc e R , rather 
than th e linear dependence upon R employed by Chiappini et al.l 
(2001). Another importa nt difference conce rns the treatment of 
star formation: in the Molla & Dfazl (I2005h model, stars form 
in two stages: first, molecular clouds condense with some effi- 
ciency out of the diffuse gas reservoir, and second, stars form 
with a second efficiency factor based upon cloud-cloud collision 
timescales . In spirit, this mimics the effect of the threshold ef- 
fect in the Chiappini et al] d2001l) model: specifically, stars may 
form only in dense regions. The relation between the star for- 
mation rate and the gas density can be approximated by a power 
law with n > 1, again, in qualitativ e agreement with the law 
employed by IChiappini et al.l d200ll) . In the halo, star forma- 
tion follows a common Schmidt-Kennicutt law with exponent 
n = 1.5. Extensive testing and tuning of the main parameters 
resulted in a grid of 440 models spanning 44 different masses 
(from dwarfs to giants, with 10 different star formation efficien- 
cies per mass model). The chemical prescriptions for SNela and 
SNell are again similar to those listed in § 2.2. 



3. Present-Day Gradients 

3.1. Radial Gradients 

In this section, the present-day radial abundance gradients of 
the MUGS and RaDES simulations are presented. We focus 
here on one MUGS (g 15784) and one RaDES galaxy (Apollo), 
which have been chosen as fiducial representatives of these two 
suites of simulations. Observational constraints on the abun- 
dance gra dient of z=0 late-ty pe galaxies may be found in, for 
example, IZaritskv et al.l d 19941) who measured a mean gradient 
of -0 .058 dex/kpc for local spiral galaxies and Ivan Zee et al.l 
d!998l) . who found a co mparable mean grad ient from their sam- 
ple (-0.053 dex/kpc). In lKewlev et al.l (2010) close galaxy pairs 
were found to have systematically shallower gradients (typically, 
-0.021 dex/kpc). In each of these cases, the gradients are in- 
ferred from gas-phase nebular emission, which provides a "snap- 
shot" of the present-day gradient, similar to that inferred from, 
for example, B-stars (i.e., stars with ages <100 Mvrs)Fl 

We employed a strict kinematic decomposition of spheroid 
and disk stars for each of the 25 simulation^! following the 
Abadi et alJ J2003) formalism. Additional (conservative) spatial 
cuts were employed to eliminate any satellite interlopers that 
might pass the initial kinematic decomposition. We define three 
age bins: young (stars born in the last 100 Myrs, to correspond 
roughly with B-stars), intermediate (stars formed 6-7 Gyr ago), 
and old (stars olders than 10 Gyr). 



10 Loose group galaxies in the RaDES suite exhibit the same qualita- 
tive flattening of metallicity gradients when compared with their 'field' 
equivalents, however the order of this difference is significantly smaller 
(<0.005 dex/kpc) than the systematic differences found between the 
RaDES and MUGS galaxies (-0.05-0.2 dex/kpc). A comprehensive 
analysis of the (subtle) systematic differences between the field and 
loose group galaxies within RaDES is forthcoming (Few et al., in prep), 
but not pursued here, simply because this difference is negligible to the 
scope of the present analysis. 

11 The kinematic decomposition e mployed for the MU GS galaxies 
differs from that used in the original Stins on et alj (2010) analysis, in 
that J-JJdrc for each star was derived self-consistently taking into ac- 
count the shape of the potential, rather than assuming spherical symme- 
try and using the enclosed mass at a given star particle's position. 
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Fig. 1. Stellar radial [Z] gradients, for disk stars in three differ- 
ent stellar populations: young (blue) is defined as stars formed 
in the last 100 Myrs, intermediate (yellow) is defined as stars 
formed 6 to 7 Gyr ago, and old (red) is defined as stars older 
than 10 Gyr. Fits to the disk are overdrawn in black; the length 
of the black line corresponds to the region of the disk used in 
the fitting (see text for details). For clarity, only two galaxies are 
shown, one from MUGS (gl5784, upper panel) and one from 
RaDES (Apollo, lower panel). 



Observational studies of radial gradients typically show 
higher m etallicities in t he inner disk relative to the outer disk 
(e.g. lRupke et a l. 2010b). As noted above, observations of exter- 
nal systems typically make use of gas-phase oxygen abundances, 
as measured from HII regions, but consistency exists between 
that tracer and others known to provide a "snapshot" of the gra- 
dient (e.g., planetary nebulae and short-lived main sequence B- 
stars). Our gas-phase and young (B-star) gradients are identical 
in amplitude and gradient, and hence in what follows, we employ 
"young stars" (those formed in the previous 100 Myr period) to 
determine the abundance gradients. 

The current RaDES sample only tracks global metallicity Z, 
but as oxygen consistently accounts for ~50% of Z, we use Z 
as a first-order proxy for oxygen, when making comparisions 
with observationsQThe version of gasoline employed for these 
MUGS runs track both O and Fe (from SNell and SNela), and 
assume ZsO+Fe; as noted earlier, this latter assumption leads to 
an ~0.2 dex underestimate of the global metallicity in the MUGS 
sample. This does not impact upon our gradient analysis, but 
does serve to explain why the RaDES and MUGS galaxies are 
offset by ~0.2 dex from one another in [Z] in the figures pre- 
sented here. 

Figure Q] shows the mass-weighted radial gradients at z-Q in 
[Z] for one MUGS galaxy (g 15784, top panel) and one RaDES 
galaxy (Apollo, lower panel). The radial gradients are calculated 
using linear fits over the noted disk regions (overdrawn in black). 



These are chosen to exclude the central region, avoiding any 
residual co-rotating bulge stars that escaped the kinematic de- 
composition. The outer edge of the disk is taken as the point at 
which the surface brightness profile of the young stars (effec- 
tively, the cold gas) deviates from an exponential. To ensure that 
an appropriate region is considered here, we have been conser- 
vative in choosing the "disk region". The gradient is robust to 
the choice of outer radius; reducing the choice of inner radius 
from 5 kpc to 2 kpc has only a +0.007 dex/kpc impact on the in- 
ferred formal gradient - i.e., the differences in gradients between 
young, intermediate, and old populations are not significantly af- 
fected. Throughout this paper we use the lAsplund et al.1 (2009) 
values for the solar metallicity. 

As one considers progressively older stellar populations (at 
the present-day), Figure Q] shows that the measured radial metal- 
licity gradient becomes progressively flatter. Such behaviour 
is not unexpected in cosmological simulations which include 
gas infall, radial flows, high velocity dispersion gas, kinemat- 
ically hot disks, and dynamical mixin g/radial migration which 
is mo re pronounced for older stars (e.g. Sanchez-Blazauez et al.1 
[20091: ?;|Pflkingion & Gibsonll2012l) . The timescale of the mix- 
ing that flattens the gradients in the MUGS and RaDES simu- 
lations is shorter than the difference between intermediate and 
old populations of stars, as evidenced by radial gradients for 
the two populations, regardless of simulation suite, being quite 
similar. The degree of flattening of the stellar abundance gra- 
dients is such that by the present day, within the simulations, 
the older stellar tracers show a flatter abundance gradient than 
the younger tracers (recall Fig 1 , re-iterating results s hown 
byJSanchez-Blazquez et al. (20091). iRahimi et all d2011l) , and 
Pilkington & Gibson] d2012h V This is counter to what is observed 
in the Milky Way when inferring gradients using y ounger plan- 
etary nebulae versus older planetary nebulae (e.g. iMaciel et al.1 
2003), but again, this is fully expected given the degree of kine- 
matic (stellar) heating within these cosmological simulations, 
and does not impact on the use of gas-phase and young-star 
probes of the gradients (both possess the expected steeper abun- 
dance gradients at early-times). Indeed, future work in this area 
can, and should, make use of this powerful constraint on migra- 
tion/heating: specifically, the fact that (empirically) older stellar 
probes today have a steeper abundance gradients than younger 
stellar probes, while extant, kinematically hot, simulations, show 
the opposite trend. 

For completeness, in Table Q] we list the present-day mass- 
weighted stellar radial metallicity gradients (d[Z]/d/?, in units 
of dex/kpc) for each of the 25 simulations employed here (col- 
umn 8). The similarity of the gradients is readily apparent, save 
for the MUGS galaxy g24334, which was included in the sam- 
ple despite its stellar fraction being dominated by accreted stars, 
rather than in s/fMStar formation (discussed further in § 3]). Its 
relatively small disk scalelength (1.0 kpc) also made fitting its 
gradient more challenging than the other MUGS disks. 



12 We have recently completed the implementation of full chemical 
evolution, including SNell, SNela, and AGB stars, within ramses - Few, 
Gibson & Courty (201 1, in prep). 



Following lSanchez^ Blazauez et alJ (1201 ll) . we examined the 
effect of applying a different weighting scheme in determining 
the mean metallicities. When examining just the young stars or 
the gas, the weighting employed has no effect upon the inferred 
gradient. However, when deriving a composite gradient making 
us e of all stars in the disk, the wei ghting can become important, 
as ISanchez-Blazquez et alJ (1201 lh suggested. We explored the 
impact of using, for example, luminosity-weighting (and log- 
weighting), by deriving the absolute magnitude of each simu- 
lated star particle, making use of its age, metallicity, and initial 
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Fig. 2. Star formation rate per unit surface area as a function of 
radius for the MUGS galaxy gl5784 (upper-left panel) and the 
RaDES galaxy Apollo (upper-right panel). We show the sim- 
ulations at four different redshifts: z=0.0, 0.5, 1.2, and 2.2, as 
noted in the inset to the upper-right panel. 1 kpc annuli are used 
along with a height cut of +5 kpc above and below the disk. 
The mass of stars formed in the last 100 Myrs is calculated for 
each annulus out to a radius of 15 kpc. The curves have been 
normalised to 1 M /Gyr/pc 2 at galactocentric radius 8 kpc. The 
bottom pa nels show the correspond ing predicted behavio ur of 
the lChiappini et alj | |2001l) (right) and lMolla & Diazl (l2005) (left) 
models. Only redshifts 0.0 and 2.2 are shown, other redshifts are 
excluded as these models evolve smoothly from z=2.2 to z=0.0. 
Two of the Molla & DH (12001 models are shown, one with 
high star formation efficiency (dashed lines) and one with low 
star formation efficiency (solid lines) . 



mass function, alongs ide thelMarigo et al.l d2008l) isoc hronesQ^I 
As expected from the ISanchez-Blazquez et alJ d201 1) analysis, 
the mean abundance shifted by ~0.1 dex depending upon the 
weighting employed, but the inferred gradient was not affected. 

The abundance gradient of young stars (or equivalently, the 
ISM) is shaped by the time evolution of the radial star forma- 
tion rate. To illustrate this we show the normalised star forma- 
tion rate per unit surface area as a function of galactocentric 
radius in Figur e [2] To match the chemical evolution models of 
Chiapmini et alJ d200ll) for the Milky Way (with the understand- 
ing that our simulations are not constructed a priori to be perfect 
replicas of the Milky Way), we normalise the star formation rate 
to have a value of 1 M /Gyr/pc 2 at a galactocentric radius of 
8kpcQ 



13 http://stev.oapd. inaf. it/cgi-bin/cmd_2 . 1 

14 For context, the 'normalised' and 'pre-normalised' star formation 
rate surface densities (at 8 kpc), for each of the simulations, are not 
dissimilar; the latter lie in the range ~1— 2 M Q /Gyr/pc 2 , save for the 
(known) discrepant MUGS galaxy g24334 (which, pre-normalised, lies 
at ~0.2 M Q /Gyr/pc 2 , reflective of the fact that its stellar content is more 
dominated by its accreted component, rather than in situ star formation. 



Each of the star formation rate profiles behave qualita- 
tively like the c l assic insid e-out chemical evolu tion models of 
IChiappini etail d2001l) and iMolla & Drazl d2005l) . in the sense 
of decreasing outwards from the inner to outer disks. An im- 
portant systematic difference between these representative sim- 
ulations is apparent though, at least at higher redshifts (l<z<2). 
Specifically, the gradient in the star formation rate per unit area 
is steeper at higher redshifts for the MUGS galaxies; it is not 
clear if this is symptomatic of a single difference between the 
MUGS and RaDES galaxies, or (more likely) a combination 
of factors including the star formation threshold, star formation 
efficiency, feedback schemes, and resolution of the respective 
simulations. Regardless, it is clear that star formation is more 
centrally-concentrated in the MUGS sample at early stages in 
the formation of the disk which unsurprisingly leads to steeper 
abundance gradients in the early disk (a point to which we return 
shortly). 

3.2. Vertical Gradients 

For completeness, as in Figs. 1 and 2, for gl5784 (MUGS) and 
Apollo (RaDES), the mass-weighted vertical stellar abundance 
gradients in the simulations are presented in Figure [3] A 'solar 
neighbourhood' is defined for each simulation as being a 2 kpc 
annulus situated at a galactocentric radius of ~2.5 disk scale- 
lengths (column 5 of Table [TJ. These radial scalelengths were 
derived from exponential fits to the stellar surface density pro- 
files. 

Classic wor k from, for exampl e, iMarsakov & Borkoval 
d2005Ll2006h and lSoubiran et alj 12008). and soon-to-be-released 
work using SDSS-SEGUE and RAVE datasets, show that verti- 
cal metallicity profiles can provide extremely effective tools for 
separating the thin disk from the thick disk. With ~300-500 pc 
softening/grid cells, we do not resolve the thin-thick disk tran- 
sition. Figure [3] shows the vertical gradient for the MUGS 
galaxy g 15784 (orange) and the RaDES galaxy Apollo (pur- 
ple), along with observatio nal d ata for the Milky Way from 
IMarsakov & Borkoval d2005l) and IMarsakov & Borkoval d2006). 
The two vertical lines show the respective resolutions of the 
MUGS and RaDES simulations. 

The vertical metallicity gradients (in their respective 'solar 
neighbourhoods') for the 25 simulations analysed here are listed 
in column 7 of TableQ] We find little variation between the simu- 
lations in question, with the typical vertical gradient lying in the 
-0.05±0.03 dex/kpc range. Only Eos, Castor, and Krios lie out- 
side this range, possessing somewhat steeper vertical abundance 
gradients. These three undergo the most extended late-time pe- 
riod of 'quiescent' evolution, as commented upon in Few et al., 
in prep. 

At face value, the vertical gradients in 

[a/HEl and [Fe/H] 

inferred from the simulations are consistent with the ob- 
served values seen in the thick disk of the Milky Way 
( — 0.05 - ~-0.08 dex/kpc). The vertical gradients in the 
Milky Way's thin disk, though, are consistently much steeper 
(where many authors find the thin dis k gradient to be be- 
tween — 0.25 -0.35 dex/kpc (e.g. ISoubiran et alj F2 008: 

Marsakov & B orkoval 120061: [Bartasiut e et all 120031: IChen et al.l 



2003)) than the results we obtain from our simulations. Our 
spatial 'resolutions' range from ~300-500 pc, and the results 



15 Here, total metallicity is used as a proxy for a in the RaDES suite, 
while oxygen is used for the MUGS and GCD+ suites; magnesium is 
used in the observational datasets described by Mar sakov & Borkoval 
l200j 120061) . 
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appear compromised on vertical scales up to ~2-3 resolution 
'elements' - i.e., any putative 'thin' disk would be (not surpris- 
ingly) unresolved. In a chemical sense, these disks are too 'hot', 
in much the same way that their ISM and stell ar populations are 
also kinematically hot (e.g. lHouse et al.ll201 ll) . 

On this issue of 'resolution', the global star formation rates 
report ed are comparative ly well converged as a function of reso- 
lution dStinson et al.l2006l §5.2.4) The most notable change with 
increasing resolution is the addition of higher redshift popula- 
tions, containing comparatively little mass, as earlier generations 
of halos are resolved. This is at least partially a result of star for- 
mation models largely being constrained to reproduce observed 
star formation rates. 

The dependence of gradients on resolution though is far less 
predictable. At our current resolution we resolve sufficient sub- 
structure and disc dynamics to capture the salient physical mech- 
anisms involved in migration. However, increasing resolution 
does resolve the physics behind migration processes better, but it 
also makes the diffusion model more localized. Equally impor- 
tantly, it is not clear to what extent the numerous processes in- 
volved in migration will interact with one another as resolution 
is increased. Taking the alternative approach of lowering reso- 
lution makes processes less likely to be captured (particularly 
substructure-induced migration), so it is not clear that conver- 
gence happens in a simple fashion. Ultimately, a definitive an- 
swer on the impact of resolution on migration requires far higher 
resolution than we are currently able to achieve and future work 
is required to address this issue. 

4. Evolution of the Radial Gradients 

While there exist a hand ful of studies of radial abundance gra- 
dients at high r edshift dJones et alj 1201 Ot ICresci et all l2010t 
I Yuan et a"Dl201 ll) . the difficulties in obtaining high resolution 
data for likely Milky Way-like progenitors has meant that the- 
oreticians have had very few constraints on their models; as 
noted earlier, inside-out galactic chemical evolution models can 
be constructed which recover the present-day gradients seen in 
the Milky Way, but they can take very different paths to get there. 
Some such models predict a steepenin g with time starting fro m 
initially inverted or flat gradients (e.g.. lChia ppini et al. (2001)), 
whil e others predict an i nitially negative gradient that flattens 
(e.g. lMolla & Dfa3 (l20Q5h \ 

To make progress in this area, we now analyse the time evo- 
lution of the gradients within our 25 simulations, supplemented 
with two classical chemical evolution models, making fits ra- 
dially at each timestep for which a clear disk could be identi- 
fied. As the disk is continually growing and evolving, we exam- 
ined each timestep visually, identifying the outer 'edge' using 
the cold gas and young stars as a demarcation point. It should 
be noted here that the kinematic decomposition used to identify 
'disk stars' in § 3. 1 and § 3.2 was not used for this component of 
our analysis. By working only with very young stars at 2-3 disk 
scalelengths, when fitting gradients at each timestep, kinematic 
decomposition of disk vs spheroid stars becomes unnecessary. 
Radial gradients were then derived by fitting typically from the 
outer edge of the disk to the inner part of the disk, where the 
inner point corresponds to the point at which the surface den- 
sity profile deviates from an exponential. Again, as we are only 
using the stars formed in the previous 100 Myrs (B-stars) at a 
given timestep, the relevant disk (rather than star-forming bulge) 
regime is not difficult to identify. 

In Figure|4] we show the time evolution of the radial gradient 
for our two 'fiducial' simulations: MUGS (gl5784, upper panel) 



g15784 [0/H] gradient = -0-06 
Apollo [Z] gradient = -0.08 
Marsakov thick disk [Mg/H] = -0.07 
Marsakov thin disk [Mg/H] = -0.16 
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g15784 [Fe/H] gradient = -0.07 
Marsakov thick disk [Fe/H] = -0.13 
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Fig. 3. The vertical gradients of disk stars in the simulations. The 
top panel shows the [Z] vertical gradient of Apollo (purple, grad 
= -0.08) with the [O/H] gradien t of g 15784 (orange, grad = 
-0.06 ) and observational data from Marsakov & Borkova ( 2005, 
2006) of [Mg/H] gradients in the thin (blue, grad = -0.16) and 
thick (green, grad = -0.07) disk of the Milky Way. The lower 
panel shows the [Fe/H] gradients of the Marsa kov & Borkoval 
d2005ll2006l) thin (grad = -0.29) and thick (grad =-0.13) disk 
data along with the g 15784 (grad = -0.07) [Fe/H] gradient. 
Overplotted vertically are the softening length of the MUGS (or- 
ange) and the minimum grid size of the RaDES (purple) simu- 
lations. The bold red lines show the region used to calculate the 
gradient. 



and RaDES (Apollo, lower panel). The gradients measured at 
each timestep are noted in the inset to each panel. Much steeper 
abundance gradients at high-redshift (z>l) are seen within the 
MUGS galaxy. Further, the offset in mean metallicity between 
the two, as already alluded to, can be traced to the manner in 
which chemistry was included in the version of gasoline em- 
ployed (i.e., the assumption that ZsO+Fe, which affects the 
mean metallicity, but not the gradient). 

In Figure |5j we show the time evolution of the [Z] gradients 
for the 4 MUGS galaxies, the gcd+ galaxy (Gall), the grape- 
SPH galaxy (KN11), and the 19 RaDES galaxies. Importantly, 
we have also derived the time evolution of the predicted gradi- 
ents for the chemical evolution models of Chiappini et al. (2001) 
and two of the Milky Way-lik e models of iMolla & Dfazl d2005l) : 
with the lMolla & Dfazl (|2005) data, the fits to determine the gra- 
dient at each timestep evolved as they did in the hydrodynam- 
ical simulations. As the disk grew, the fits were made at larger 
radii, to exclude the central region. From the earliest timestep to 
the latest the fitted region shifts ~3 kpc in radius (reflecting the 
growth of the disk over the timescales under consideration). The 
Chiappini et alj d2001l) data were fit over the radial range 4 to 
8 kpc at each timestep, reflecting the fewer relevant an nuli avail- 
able over which to make the fit. Chia ppini et al.l (1200 ll) fit their 
gradients to the same chemical evolution models over a broader 
radial range (4-14 kpc), but our interests here are restricted to 
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the inner disks of these models, where the star formation density 
threshold is less important i n shapi ng the metallicity gradient. 

For the iMolla & Dfazl (120051) models, we show a low- 
efficiency (28,8) and high-efficiency (28,2) example, (where 
model 28 corresponds to a circular velocity of ~200 km/s and 
the efficiency factors correspond to the combined efficiency 
of molecular cl o ud for mation and cloud-cloud collisions). The 
Chiappini et al.1 (1200 lh and , to a lesser extent, the high effi- 
ciency Molla & Dfazl d20 05) models (at least since z~l) steepen 
with time Conversely, the RaDES sample (represented by the 
purple hatched region, which encompasses l<x of the gradient 
values at a given redshift) shows a mild flattening with time, 
more in keeping wi th full time evolution of the high efficiency 
IMolla & Dfazl (120051) model. The MUGS sample shows not only 
steeper gradients as a whole at z>l (except for g24334, to which 
we return below), but also three of the four show the more sig- 
nificant degree of flattening alluded to in relation to Fig |4l this 
degree of flattening is more dramatic than that seen in any of the 
RaDES galaxies or the che mical evolution mode ls (except for 
the low efficiency models of Moll a & Dfazl (l2005l) )Fl 

Shown also in Fig |5] are the typica l gradients encountered in 
nearby isolated (jZa ritskv et al. ( 1994); blue asterisk) and inter- 
acting dKgwlexeray (1201(3); red asterisk) disk galaxies (offset at 
z=0, for clarity, in Fig[5]l. The black asterisk at redshift z~ 1 .5 cor- 
responds to the recent determination of a stee p metallicity gradi - 
ent in a high-redshift grand design spiral by I Yuan et al.l (1201 lh . 
While intriguing, it is important to bear in mind that one should 
not necessarily make a causal link between these disparate data 
points; until a statistical sa mple of high-redsh ift gradients has 
been constructed, linking the lYuan et al.l (1201 lh point with those 
at low-redshift should be done with caution. 

For this latter reason, we have also included one MUGS 
galaxy (g24334) in our analysis (red curve: Fig [5j that does not 
have a present-day gradient consistent with the typical late-type 
spiral. We chose to include it, in order for the reader to see one 
example of a disk which possesses a steep gas-phase abundance 
gradie nt at high-redshift, comparable in slope to the Yuan et al. 
(2011) observation, but one which does not evolve in time to 
resemble the shallower slopes seen in nature today. g24334 dif- 
fers from the other MUGS galaxies, in the sense that the fraction 
of its stellar population born 'in situ', as opposed to 'accreted', 
is significantly lower. Further, its disk is less extended than the 

16 The Chiappini et al. ( 2001) models have gradients which are mildly 
inverted at high-redshift ( — 1-0.02 dex/kpc at redshift z~2); t his works in 
the sa me direction as the inverted gradients observed by Cres ci et al] 
d2010h at z~3, albeit the gradients claimed by the latter are signifi- 
cantly more inverted (i.e., — 1-0. 1 dex/kpc) than encountered in any of 
the simulations or chemical evolution models. It is important to remem- 
ber though that the AMAZE/LSD samples at z~3.3 are (a) primarily 
Lyman-Break Galaxies with star formation rates (~ 100-300 M G /yr) 
well in excess of that expected for Milky Way-like progenitors, and 
are not likely ideal progenitors against which to compare these simu- 
lations or chemical evolution models, and (b) in none of the current 
simulations are we able to unequivocally identify stable rotationally- 
supported disks, like those compiled by AMAZE/LSD. We require tar- 
geted simulations with much higher resolution at high-redshift than we 
have access to here, and tuned to be more representative of high-redshift 
Lyman-break galaxies, before commenting further on this potentially 
interesting constraint. 

17 It is worth noting that no obvious trend is seen when comparing 
the field and group galaxies in the RaDES sample. This is perhaps 
attributable to our selection criteria; by removing strongly interacting 
galaxies (at or near a pe ricentre passage), the s ort of systematic d iffer- 
ences seen in the work of Rupke et al. (2010a b); Perez et all d201 ll) . for 
example, would not be encountered here. 



other Milky Way-analogs and its abundance gradient was de- 
rived at ~0.5x disc scalelengths, where the grad ient is more ro- 
bust to interaction-induced flattening (e.g. lPerez et al.ll201 lh . 

These differences are ultimately traced to the underlying 
treatment of star formation and feedback within the simula- 
tions; for example, the MUGS galaxies have a higher star for- 
mation threshold than the RaDES suite (1 cm" 3 vs 0.1 cm 3 ). 
As su ch, both the MUGS s ample and the low efficiency mod- 
els of Molla & Dfazl (1200 5) preferentially form stars in the in- 
ner disk where the densities are higher; the RaDES galaxies and 
the remaining chemical evolution models, with the lower thresh- 
old, have star formation occurring more uniformly throughout 
the early disk. Further, both MUGS and RaDES employ a stan- 
dard blast-wave for malism for energy deposition into the ISM 
(IStinson et al.1 [2006). but the latter imposes a minimum blast 
wave radius of 2 grid cells, which means that ejecta is in some 
sense more "localised" in the MUGS simulations (for the same 
SN energy, the RaDES blast waves are ~2-3x larger); distribut- 
ing energy (and metals) on larger radial scales can result in a 
more uniform (i.e., flattened) metallicity distribution. The trend 
of Gall lies somewhat between the extremes of MUGS and 
RaDES, which can be traced to the fact that Gall uses a lower 
star formation threshold density (0.1 cm" 3 ), and almost negli- 
gible feedback, resulting in more localised metal enrichment. 
KN11 also lies very close to the MUGS fiducial (g 15784) in 
terms of the temporal evolution of its abundance gradient; both 
employ high SNe feedback efficiencies, albeit on different spa- 
tial scales (a density-dependent blast wave radius in the case 
of gl5784 and a fixed 1 kpc radius in the case of KN11) and 
with different star formation prescriptions (a 1 cm" 3 star forma- 
tion density threshold in the case of gl5784 and an absence of a 
threshold for KN11). Note that although these hydrodynamical 
simulations experience different merger histories, the metallic- 
ity gradients are more affected by the recipe of sub-grid physics. 
This is highlighted by our large samples of simulations gener- 
ated with different codes. 

As detailed in § 2.5, Chiappini et al] (1200 lh use a two infall 
model; at early times the infall of primordial gas is rapid and 
independent of galactocentric radius, while at later times, gas is 
assumed to fall preferentially on the outer regions of the disk, 
causing a steepening of the gradient with time. The radial de- 
pendence of this disk infall timescale is fairly gentle ( linear with 
increasing radius); on the other hand, Moll a & Dfazl d2005h cal- 
culate the overall infall rate as a function of the mass distribution 
and rotation of the galaxy, and assume a much stronger radial 
dependence for the infall timescale. Specifically, the inner disk's 
infall timescale is much more rapid than that of Chiappini et al.l 
(2001), while the outer disk's infall timescale is much longer. In 
combination, the gradient tends to flatten with time (particularly 
for their low efficiency models). 

We find clear evidence of inside-out formation in the star for- 
mation profiles at different redshifts. Starting from an initially 
concentrated distribution, this flattens with time to the present- 
day, where star formation is more extended (and close to con- 
stant) over a large fraction of the disk (Fig|2|. The radial depen- 
dence of star formation rate to infall rate sets the magnitude of 
the abundance gradient dChiappini et al.ll2001l) ; a stronger radial 
dependence resulting in a steeper gradient. Such a configuration 
appears to come about naturally in the MUGS simulations, due 
in part to their higher star formation rate density threshold and 
perhaps the higher star formation efficiency and more localised 
chemical/energetic feedback. This contributes to the steeper gra- 
dients seen at early times in these simulations, relative to the 
other models. The RaDES galaxies behave more like the high 
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Fig. 5. The derived radial [Z] gradient as a function of redshift. 
Here, we have used 1 1 different redshifts and measured the radial 
gradient of the young stars (stars formed in the last 100 Myrs at 
each step) in the disk at that time. We examined the disks at each 
redshift, to determine the appropriate galactocentric radius over 
which to measure the gradients (see text for details). Four MUGS 
galaxies (gl5784 (orange diamonds); g24334 (red diamonds); 
g422 (black diamonds); gl536 (gr een diamonds)) are shown, 
along with Gall (blue squares) from Rahimi et alj ( 201 H). K N1 1 
(cyan plus symbols) from Kobayashi & Nakasato (2011), and 
the 19 RaDES galaxies (denoted by the purple hatched area 
showing the region encapsulating 1 cr of the gradients measured 
at a given redshift). The two chemical evolution models are over- 
laid for completeness: Chiappini (black dot dashed crosses), and 
Molla high efficiency (black dashed triangles) and low effiency 
(black dotted triangles). The black asterisk correspo nds to the 
result from one lensed grand design spiral at z~1.5 dYuan et alj 
1201 ll) . the blue asterisk to th e typical gradient inferred in nearby 
spirals (IZaritskv et al.lll994h . and the the red asterisk to th e typi- 
cal gradient seen in interacting disks dKewlev et alj|2010l) ; these 
latter local points are offset slightly at z=0, for clarity. 



radius [kpc] 



Fig. 4. The radial [Z] gradients of young stars in gl5784 (top 
panel) and Apollo (bottom panel). The different colors corre- 
spond to different redshifts running from z=0 (black) to z=2.2 
(orange), illustrating the time evolution of the abundance gradi- 
ents in both simulations. Note the more dramatic flattening of the 
MUGS (gl5784) relative to that of RaDES (Apollo). The fitted 
gradients were not done in an 'automated' fashion; we examined 
each timestep's surface density, kinematic, and abundance pro- 
files, to take into account the growth of the disk and identify the 
'cleanest' disk region within which to determine the gradient. 

efficiency model of Molla & Dfazl d2005l) . It should be noted 
however that despite the significant differences seen in the early 
stages of these galaxies' evolution, the star formation distribu- 
tion in the majority of these simulations is very similar at the 
present day. 



5. Summary 

This work provides evidence in support of the imposed inside- 
out disk growth paradigm adopted within chemical evolution 
models; this growth is a natural outcome of both Eulerian and 
Lagrangian hydrodynamical simulations of disk galaxy forma- 
tion within a cosmological context. We have examined how this 
inside-out growth impacts on the magnitude and evolution of 
abundance gradients in these galaxies, using a suite of simula- 
tions and models which were calibrated to recover the present- 
day shallow gradients observed in late-type spirals. This is not 
meant to be a comprehens ive, systematic, examination of sub- 
grid physics, in the vein of Wiers ma et alj (1201 ll) . for example; 
instead, we have taken (in some sense) the 'best' Milky Way -like 
simulations from several groups, using different codes, differ- 
ent initial conditions, and different assembly histories, and con- 
ducted a 'blind' experiment on the outputs, to quantify how the 
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gradients evolved to the imposed boundary condition of a shal- 
low present-day gradient. Our findings include the following: 

1. All galaxy models and simulations described in this work 
exhibit inside-out formation of the disk with varying de- 
grees of centrally-concentrated star formation at early times 
(Figure |2J. The evolving radial star formation rate depen- 
dence directly influences the resulting metallicity gradient; 
put another way, the signature of the star formation profile 
is embedded within the gradient of the young stars at each 
timestep. This signature though is diluted on the timescale 
of a few Gyrs. This is reflected in the differing gradients 
at the present-day between old and young stars (Figure Q]); 
young stars at high-redshift within the MUGS sample (and 
observationally, it would appear, tentatively) form with a 
steep metallicity gradient, while those same stars today (now, 
old) have a fairly flat metallicity gradient (see Pilkington & 
Gibson 2012). 

2. Within the suite of 25 cosmological hydrodynamical simula- 
tions the derived vertical abundance gradients are compara- 
ble to those observed locally in the Milky Way's thick disk. 
The resolution is, however, not sufficient to discriminate be- 
tween thin and thick disks. 

3. The evolution of simulated metallicity gradients depends 
strongly on the choice of sub-grid physics employed and 
as such the magnitude and direction of its evolution de- 
pends critically upon the specific details of the recipes im- 
plemented. While it is difficult to disentangle the behaviour 
of the star formation profile a priori, it is clear that simu- 
lated galaxies with more centrally-concentrated star forma- 
tion have initially steeper abundance gradients. These are 
more consistent with the (albeit limited) observat ion of high 
redsh ift normal Grand Design spiral galaxies dYuan et al.l 
12011b . 

4. All the models and simulations tend to similar present-day 
abundance gradients, despite the diversity at earlier times, 
save for g24334 (which was chosen specifically in violation 
of the imposed shallow present-day gradient boundary con- 
dition, for illustrative purposes). In almost every case this 
requires the gradient to flatten with time, the exception be - 
ing the chemical evolution model of Chiapp ini et al.l (120011) . 
This model starts with an initially positive gradient that is 
independent of its halo phase. The gradient then inverts to 
become negative, with a gradient similar to other chemical 
evolution models. 

5. The diversity of the evolution of metallicity gradient is for 
the first time highlighted by our large sample of both hydro- 
dynamical simulations and chemical evolution models. Our 
results indicate that observations of the metallicity gradient 
for disk galaxies at different redshifts and that for the dif- 
ferent age populations in the Galaxy are key to reveal the 
formation processes of disk galaxies and better constrain the 
sub-grid physics implemented with all the codes sampled. 

Future work in this area will see us employ a finer temporal ca- 
dence, in order to better track the precise influence of merger 
events on the abundance gradients (both the magnitude of the 
effect and the timescale for re-establishing a stable abundance 
gradient). This study will also yield a deeper understanding of 
how the non-linear processes of star formation and feedback in- 
fluence systematic differences between the various simulations 
presented here. We are near completion of a major upgrade 
to ramses which will allow us to re-simulate the RaDES suite 
with a broad spectrum of chemical elements, including those 



from SNell, SNela, and AGB stars. With ongoing and future 
large scale spectroscopic surveys and missions such as RAVE, 
APOGEE, SEGUE, HERMES, LAMOST, and Gaia, providing 
detailed information on the phase and chemical space signatures 
of the Milky Way and beyond, such a chemodynamical explo- 
ration will be both timely and critical for understanding the ori- 
gin and evolution of abundances in galaxies, and their link to the 
underlying physics of galaxy formation. 
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